function  x = qmrirecon( val, err, varargin )
% create regularised reconstruction of quantitative MRI map 
% by minimizing F(x) = (f-x)' B^-1 (f-x) + alpha* xL'Lx + beta x'x
% with: 
% f = the estimated values (the T-map)
% x = the regularised values
% B = diagonal matrix with the error values (the errormap) 
% alpha = weighting factor for laplacian smoothing
% beta = weighting factor for high value removal
% solution: dF/dx = 0 -> -2 B^-1 (f-x) + 2 alpha L'Lx + 2 beta x = 0
% -> (B^-1 + alpha L'L + beta) x = B^-1 f

% Usage:
% x = qmrirecon( val, err, alpha, beta, err_is_stdev )
%
% where:
% val = the estimated value (dcm filename, or matrix)
% err = the cramerrao bound (standard deviations) (dcm filename, or matrix)
% alpha = weighting factor for laplacian; default 0.001 if omitted
% beta = weighting factor for identity matrix; default 0.001 if omitted
% err_is_stdev = the provided err map is actually the sqrt of the
% cramer-rao bound; default true if omitted.
% x = the regularised reconstruction of the map.
%
% Example usage:
% first set the directory to where your images are
%    cd('C:\User\MRIdata')
% then execute
%    qmrirecon('0006T1map.dcm','0006T1errormap.dcm', 0.01,0.0001)

%
% author: Stefan Klein, 14-02-2011.
% modified by Henk Smit, questions: h.smit@erasmusmc.nl

name = val(1:(size(val,2)-4));
if ( ischar(val) )
  val = double( dicomread( val ) );
end
if ( ischar(err) )
  err = double( dicomread( err ) );  
end

if numel(varargin) > 0
  alpha = varargin{1};
else
  alpha = 0.001;
end

if numel(varargin) > 1
  beta = varargin{2};
else
  beta = 0.001;
end
   
if numel(varargin) > 2
  err_is_stdev = varargin{3};
else
  err_is_stdev = true;
end

if err_is_stdev
    err = err.^2;
end


lapfilt=[0 1 0; 1 -4 1; 0 1 0];
sizex=size(val);
eps=1e-16;
eps=0.1;
n = numel(val);

inverr = 1./(err(:)+eps);

%x=pcg( @(x)Ax(x, inverr, lapfilt, alpha, sizex) , r2(:).*inverr, 1e-6, 100, spdiags(err(:)+1, 0,n,n) );
x=pcg( @(x)Ax(x, inverr, lapfilt, alpha, beta, sizex) , val(:).*inverr, 1e-6, 500 );

x=reshape( x, sizex );

%clim=[min(x(:)) max(x(:))];
clim=[0 20];

figure;
imagesc(val);
axis image;
set(gca,'CLim', clim);
colormap('gray');
colorbar;
title( 'value' );
set(gca,'XTick', []);
set(gca,'XTickLabel', []);
set(gca,'YTick', []);
set(gca,'YTickLabel', []);

figure;
imagesc(err);
axis image;
set(gca,'CLim', clim);
colormap('gray');
colorbar;
title( 'error' );
set(gca,'XTick', []);
set(gca,'XTickLabel', []);
set(gca,'YTick', []);
set(gca,'YTickLabel', []);

figure;
imagesc(x);
axis image;
set(gca,'CLim', clim);
colormap('gray');
colorbar;
title( 'regularised' );
set(gca,'XTick', []);
set(gca,'XTickLabel', []);
set(gca,'YTick', []);
set(gca,'YTickLabel', []);

dicomwrite(int16(x),[name,'_Regularized','.dcm']);
dicomwrite(int16(val),[name,'_orig','.dcm']);



function y = Ax(x, inverr, lapfilt, alpha, beta, sizex)

y= conv2( reshape(x, sizex), lapfilt, 'same' );
y= alpha*conv2( y, lapfilt, 'same' );
y=y(:)+x.*(inverr+beta); 

